Monte Carlo Simulations and LLN

Monte Carlo Simulations

General Idea of Monte Carlo Analysis (for this course)

• Use computer simulation to stochastically approximate some expected value.
  • Since expected value of an indicator variable is a probability, same method includes stochastically approximating probabilty of some event.
  • Relies critically on the Law of Large Numbers.

Monte Carlo Simulations

General Idea of Monte Carlo Analysis (for this course)

• Use computer simulation to stochastically approximate some expected value.
  • Use when deriving an exact solution is difficult or infeasible.
  • Often not computationally efficient.

Monte Carlo Simulations

MC simulations used to:

• Simulate expected payoff from some strategy;
• Study properties of estimators and inference procedures;
• Study quality of asymptotic approximations . . .

Monte Carlo Simulations

MC simulations used to:

• Implement some inference procedures
  • Have already done so with randomization inference,
  • will return to with bootstrap infernece;
• Implement some estimation procedures
  • for example, creation of training and test samples in machine learning.
  • we will use later in this course for K-fold cross-validation.

pass line bet in craps

• Otherwise, shooter keeps rolling, until roll same sum as original sum (shooter wins) or sum is 7 (casino wins).

What is the probability that shooter wins?

pass line bet in craps

• In this example, could use probability theory to derive that the probability of winning is \(244/495 \approx 0.4929\).

  • See, e.g., here.

  • But might be easier to use MC simulation than to do the derivation.

• Not feasible to derive solution in many examples.

Win on first roll

• Let’s start with simpler problem:
  what is the probability that shooter wins on first roll?

  • Can derive that probability equals \(2/9 \approx 0.2222\)

  • Now we use Monte Carlo simulation to approximate that probability.

Rolling a Die

• First step: simulating one roll of a die:

  • R code: sample(1:6,1) which randomly drawing an element from {1,2,3,4,5,6}, each element equally likely, one time.

sample(1:6,1)

[1] 6

Rolling a Die

• Rolls are random, each time can get different results:

sample(1:6,1)

[1] 5

sample(1:6,1)

[1] 3

sample(1:6,1)

[1] 6

sample(1:6,1)

[1] 6

Are results really random?

• Computers are deterministic.

  • To create truly random numbers, can use a physical process, e.g., from nuclear decay or radioisotopes.

  • Not necessary for our purposes, and not what R is doing…

  • R is generating pseudo-random numbers

R generating pseudo-random numbers

• R is generating pseudo-random numbers

  • starting point is a seed, in R determined by default by computer’s real time clock (time to the nanosecond that you run process), along with process ID number;
  • algorithm applied to seed to obtain pseudo-random numbers.
  • not actually random, as deterministic given seed, but close enough to random for our purposes.

Replication?

• You may obtain different pseudo-random numbers each time you run your code, neither you nor others can replicate your results.

• For replication purposes, can set.seed

set.seed(123)
sample(1:6,1)
sample(1:6,1)
sample(1:6,1)

[1] 3

[1] 6

[1] 3

set.seed(123)
sample(1:6,1)
sample(1:6,1)
sample(1:6,1)

[1] 3

[1] 6

[1] 3

Rolling Two Dice

• Now simulating rolling two dice:

  • sample(1:6,2,replace=TRUE) randomly draws one element from {1,2,3,4,5,6}, each element equally likely, with replacement, two times.

set.seed(123)
sample(1:6,2,replace=TRUE)

[1] 3 6

set.seed(123)
sum(sample(1:6,2,replace=TRUE))

[1] 9

Approximate Prob. Win on 1st Roll

• Shooter wins on first roll if sum of dice is 7 or 11.

• How to approximate probability winning first roll?

• Simulate many rolls, find what fraction of rolls results in a sum of 7 or an 11.

Approximate Prob. Win on 1st Roll

• More formally
  • Let \(X\) denote indicator variable for winning on first roll.
  • Let \(p = \Pr[X=1]\), probability win.
  • Let \(X_1, . . . ,X_n\) denote \(n\) replications of \(X\).
  • Let \(\bar{X}_n\) sample mean, i.e., fraction of replications with \(X_i = 1\).
  • Use \(\bar{X}_n\) to approximate \(\mathbb{E}[X]=\Pr[X=1]\).

Win first roll

• One simulation of whether win on first roll:

sum(sample(1:6,2,replace=TRUE)) %in% c(7, 11) 

[1] FALSE

Approximate Prob. Win on 1st Roll

• Replicate that simulation 10 times

set.seed(123)
replicate(10,sum(sample(1:6,2,replace=TRUE)) %in% c(7, 11) )

 [1] FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE

set.seed(123)
mean(replicate(10,sum(sample(1:6,2,replace=TRUE)) %in% c(7, 11) ))

[1] 0.1

Approximate Prob. Win on 1st Roll

• True probability win on first roll: \(2/9 \approx 0.222\).

• Fraction of 10 simulations where won: 0.1.

  • Not a very good approximation.

Approximate Prob. Win on 1st Roll

• Additional problem:
  Can get very different approximation if rerun code:

mean(replicate(10,sum(sample(1:6,2,replace=TRUE)) %in% c(7, 11) ))

[1] 0.1

mean(replicate(10,sum(sample(1:6,2,replace=TRUE)) %in% c(7, 11) ))

[1] 0.4

mean(replicate(10,sum(sample(1:6,2,replace=TRUE)) %in% c(7, 11) ))

[1] 0.2

Approximate Prob. Win on 1st Roll

• Additional problem:
  Can get very different approximation if rerun code:

  • Note that \(\mathbb{E}(\bar{X}_n)=p\), but \(\mbox{Var}(\bar{X}_n)= p*(1-p)/N\), where \(p=2/9\).

Stochastic Approximation:

• Consider increasing number of replications

# mean of 100 rolls of a die
mean(replicate(100,sum(sample(1:6,2,replace=TRUE))
               %in% c(7, 11) ))

[1] 0.25

# mean of 1000 rolls of a die
mean(replicate(1000,sum(sample(1:6,2,replace=TRUE)) 
               %in% c(7, 11) ))

[1] 0.203

# mean of 10000 rolls of a die
mean(replicate(10000,sum(sample(1:6,2,replace=TRUE))
               %in% c(7, 11) ))

[1] 0.2157

# mean of 100000 rolls of a die
mean(replicate(100000,sum(sample(1:6,2,replace=TRUE)) 
               %in% c(7, 11) ))

[1] 0.22289

Stochastic Approximation:

• Plotting results for \(n=1\) to \(=1000\).

Stochastic Approximation:

• What happens if repeat this MC simulation?

Stochastic Approximation:

• What if we repeat this simulation 100 times?

Stochastic Approximation:

Use Theory to Understand?

• Recall Convergence in Probability

  • We say a sequence \(\bar{X}_n\) converges in probability to \(\mu\), written \(\bar{X}_n \stackrel{p}{\rightarrow} \mu\), if, for every \(\epsilon>0\), \[ \lim_{n \rightarrow \infty} \Pr[ | \bar{X}_n - \mu | \ge \epsilon] = 0.\]

    • intuitive meaning?

  • If \(\bar{X}_n \stackrel{p}{\rightarrow} \mu\), we say that \(\bar{X}_n\) is a consistent estimator of \(\mu\).

Law of Large Numbers

• Recall LLN for i.i.d. data.

  • Let \(X_1, . . . , X_n\) denote i.i.d. random variables, with \(\mathbb{E}[X_i]=\mu,\) finite. Then \(\bar{X}_n \stackrel{p}{\rightarrow} \mu\).

• When \(X_i\) an indicator variable, \(\mathbb{E}[X_i]=p\) with \(0 \le p \le 1\), so \(\mathbb{E}[X_i]=\mu\) always exists, finite. Thus, by LLN, \(\bar{X}_n \stackrel{p}{\rightarrow} p\), i.e., \(\bar{X}_n\) approximates \(p\) arbitrarily well for \(n\) sufficiently large.

Original Problem

• What is the probability that shooter wins pass line bet in craps?

• First write function to simulate pass line bet.

Flowchart for Craps

flowchart LR
D -->| &nbsp No &nbsp | F( &nbsp Roll <br> Again &nbsp )
subgraph  Roll Again &nbsp
F(Roll Again  &nbsp) --> G[Roll     7? &nbsp ]  
  G -->|Yes &nbsp | H{ &nbsp Lose &nbsp }
  G -->|No &nbsp | I[Roll = <br>    1st   roll? &nbsp]
  I -->|Yes &nbsp | J{ &nbsp Win &nbsp }
  I -->|No &nbsp| F
  end
  subgraph 1st Roll &nbsp
  A(&nbsp Roll <br> Dice &nbsp ) --> B[&nbsp  Roll <br>  &nbsp 7 or 11? &nbsp ]
  B -->| &nbsp Yes &nbsp | C{ &nbsp Win &nbsp }
  B -->| &nbsp No &nbsp | D[&nbsp  Roll 2 <br>  &nbsp  3   or 12? &nbsp ]
  D -->| &nbsp Yes &nbsp | E{ &nbsp Lose &nbsp }
end

Function to simulate pass line bet

f.craps <- function(){
    # simulate first roll of dice
    sum0 <- sum(sample(1:6,2,replace=TRUE))
    # determine if win/lose on first roll
    # determine outcome if don't win/lose on first roll
    return(Win)
}

Function to simulate pass line bet

f.craps <- function(){
    # simulate first roll of dice
    sum0 <- sum(sample(1:6,2,replace=TRUE))
    # determine if win/lose on first roll
    if (sum0  %in% c(7, 11)){
             Win <- 1
    } else if (sum0 %in% c(2, 3, 12)){
             Win <- 0 
    } else{
    # determine outcome if don't win/lose on first roll
    }
    return(Win)
}

Function to simulate pass line bet

f.craps <- function(){
    # simulate first roll of dice
    sum0 <- sum(sample(1:6,2,replace=TRUE))
    # determine if win/lose on first roll
    if (sum0  %in% c(7, 11)){
             Win <- 1
    } else if (sum0 %in% c(2, 3, 12)){
             Win <- 0 
    } else{
    # determine outcome if don't win/lose on first roll
    while(TRUE){
    # keep rolling until someone wins
    }
    }
    return(Win)
}

Function to simulate pass line bet

f.craps <- function(){
    # simulate first roll of dice
    sum0 <- sum(sample(1:6,2,replace=TRUE))
    # determine if win/lose on first roll
    if (sum0  %in% c(7, 11)){
            Win <- 1
    } else if (sum0 %in% c(2, 3, 12)){
            Win <- 0 
    } else{
    # determine outcome if don't win/lose on first roll
    while(TRUE){
    # keep rolling until someone wins
      sum1 <- sum(sample(1:6,2,replace=TRUE))
      if (sum1 == sum0){
            Win <- 1
            break
        } else if (sum1 == 7){
            Win <- 0
            break
        }
    }
    }
    return(Win)
}

Simulate Prob Win Pass Line Bet

sim.prob <- mean(replicate(100000,f.craps()))
sim.prob 

[1] 0.49262

• MC simulated probabilty: 0.49262
• True prob: - \(244/495 \approx 0.4929\).

General Procedure

From some function \(f\) ,

• Let \(X\) denote random variable (or random vector).

• Let \(X_1, . . . ,X_R\) denote \(R\) replications of \(X\).

• Construct \(f(X_1), . . . ,f(X_R)\).

• Let \(\overline{f(X)}_R\) sample mean of \(f(X_j)\) across \(R\) replications.

• Use \(\overline{f(X)}_R\) to approximate \(\mathbb{E}[f(X)]\).

Examples:

• Approximate \(\mathbb{E}[X]\) for \(X \sim N(0,1)\).

mean(rnorm(100000))

[1] 0.001695661

f.1 <- function(x){max(0,x)} 
mean(replicate(100000,f.1(rnorm(1))))

[1] 0.3995883

# mean of 10000 replications
f.2 <- function(x){(-1.96 <x & x<1.96)} 
mean(replicate(100000,f.2(rcauchy(1))))

[1] 0.70044

Examples:

Let \(\bar{Y}_{10}\) be the mean of \(10\) fair coin flips. Using \(100,000\) replications,

• approximate \(\Pr[ | \bar{Y}_{10} -0.5 | < 0.1]\)

mean(replicate(100000,abs(mean(rbinom(10,1,p=0.5))-0.5)<0.1))

[1] 0.65691

• approximate \(\Pr[ | \frac{1}{\bar{Y}_{10}} -2 | < 0.1]\)

f.inv <- function(x){1/x} 
mean(replicate(100000,abs(f.inv(mean(rbinom(10,1,p=0.5)))-2)<0.1))

[1] 0.2464

Examples:

Let \(\bar{Y}_{100}\) be the mean of \(100\) fair coin flips. Using \(100,000\) replications,

• approximate \(\Pr[ | \bar{Y}_{100} -0.5 | < 0.1]\)

mean(replicate(100000,abs(mean(rbinom(100,1,p=0.5))-0.5)<0.1))

[1] 0.96573

• approximate \(\Pr[ | \frac{1}{\bar{Y}_{100}} -2 | < 0.1]\)

f.inv <- function(x){1/x} 
mean(replicate(100000,abs(f.inv(mean(rbinom(100,1,p=0.5)))-2)<0.1))

[1] 0.38007

Examples:

Let \(\bar{Y}_{1000}\) be the mean of \(1,000\) fair coin flips. Using \(100,000\) replications,

• approximate \(\Pr[ | \bar{Y}_{1000} -0.5 | < 0.1]\)

mean(replicate(100000,abs(mean(rbinom(1000,1,p=0.5))-0.5)<0.1))

[1] 1

• approximate \(\Pr[ | \frac{1}{\bar{Y}_{1000}} -2 | < 0.1]\)

f.inv <- function(x){1/x} 
mean(replicate(100000,abs(f.inv(mean(rbinom(1000,1,p=0.5)))-2)<0.1))

[1] 0.88559

Examples:

Let \(\bar{Y}_{10000}\) be the mean of \(10,000\) fair coin flips. Using \(100,000\) replications,

• approximate \(\Pr[ | \bar{Y}_{10000} -0.5 | < 0.1]\)

mean(replicate(100000,abs(mean(rbinom(10000,1,p=0.5))-0.5)<0.1))

[1] 1

• approximate \(\Pr[ | \frac{1}{\bar{Y}_{10000}} -2 | < 0.1]\)

f.inv <- function(x){1/x} 
mean(replicate(100000,abs(f.inv(mean(rbinom(10000,1,p=0.5)))-2)<0.1))

[1] 1

Examples:

Let \(\bar{Y}_{10}\) be the mean of \(10\) fair coin flips. Using \(100,000\) replications,

• approximate \(\Pr\left[~ \left | \sqrt{N} \frac{\bar{Y}_{10} -0.5}{\sqrt{\bar{Y}_{10} * (1-\bar{Y}_{10})}} \right | < 1.96\right]\)

f.tstat0 <- function(n){
  phat <- mean(rbinom(n,1,p=0.5))
  sqrt(n)*(phat - 0.5)/sqrt(phat * (1-phat))
}
  
mean(replicate(100000,abs(f.tstat0(10))<1.96))

[1] 0.891

Examples:

Let \(\bar{Y}_{100}\) be the mean of \(100\) fair coin flips. Using \(100,000\) replications,

• approximate \(\Pr\left[~ \left | \sqrt{N} \frac{\bar{Y}_{100} -0.5}{\sqrt{\bar{Y}_{100} * (1-\bar{Y}_{100})}} \right | < 1.96\right]\)

f.tstat0 <- function(n){
  phat <- mean(rbinom(n,1,p=0.5))
  sqrt(n)*(phat - 0.5)/sqrt(phat * (1-phat))
}
  
mean(replicate(100000,abs(f.tstat0(100))<1.96))

[1] 0.94338

Examples:

Let \(\bar{Y}_{1000}\) be the mean of \(1000\) fair coin flips. Using \(100,000\) replications,

• approximate \(\Pr\left[~ \left | \sqrt{N} \frac{\bar{Y}_{1000} -0.5}{\sqrt{\bar{Y}_{1000} * (1-\bar{Y}_{1000})}} \right | < 1.96\right]\)

f.tstat0 <- function(n){
  phat <- mean(rbinom(n,1,p=0.5))
  sqrt(n)*(phat - 0.5)/sqrt(phat * (1-phat))
}
  
mean(replicate(100000,abs(f.tstat0(1000))<1.96))

[1] 0.94661

Examples:

Let \(\bar{Y}_{10}\) be the mean of \(10\) fair coin flips. Using \(100,000\) replications,

• approximate \(\Pr\left[~ \left | \sqrt{N} \frac{\bar{Y}_{10} -0.7}{\sqrt{\bar{Y}_{10} * (1-\bar{Y}_{10})}} \right | < 1.96\right]\)

f.tstat0 <- function(n){
  phat <- mean(rbinom(n,1,p=0.5))
  sqrt(n)*(phat - 0.7)/sqrt(phat * (1-phat))
}
  
mean(replicate(100000,abs(f.tstat0(10))<1.96))

[1] 0.81663

Examples:

Let \(\bar{Y}_{100}\) be the mean of \(100\) fair coin flips. Using \(100,000\) replications,

• approximate \(\Pr\left[~ \left | \sqrt{N} \frac{\bar{Y}_{100} -0.7}{\sqrt{\bar{Y}_{100} * (1-\bar{Y}_{100})}} \right | < 1.96\right]\)

f.tstat0 <- function(n){
  phat <- mean(rbinom(n,1,p=0.5))
  sqrt(n)*(phat - 0.7)/sqrt(phat * (1-phat))
}
  
mean(replicate(100000,abs(f.tstat0(100))<1.96))

[1] 0.01728

Examples:

Let \(\bar{Y}_{1000}\) be the mean of \(1000\) fair coin flips. Using \(100,000\) replications,

• approximate \(\Pr\left[~ \left | \sqrt{N} \frac{\bar{Y}_{1000} -0.7}{\sqrt{\bar{Y}_{1000} * (1-\bar{Y}_{1000})}} \right | < 1.96\right]\)

f.tstat0 <- function(n){
  phat <- mean(rbinom(n,1,p=0.5))
  sqrt(n)*(phat - 0.7)/sqrt(phat * (1-phat))
}
  
mean(replicate(100000,abs(f.tstat0(1000))<1.96))

[1] 0

Explanation? Interpretation

Switching to blackboard. . .